* Settings
version 16
do "$SSDIMed/scripts/_auxiliary/_project_settings.do"


* -------------------------------------------------------------------------------------------------
* Prepare data for figures
* -------------------------------------------------------------------------------------------------

* DI entrant outcomes, by x, where x = 
*   i.covstart_year          = calendar year of Medicare entry
*   c.unemprateatapp         = unemployment rate at time of application
*   i.unemp_q20              = ventiles of unemprateatapp
*   i.age_year_covstart_fill = age in years at Medicare entry

* Sample
*   main: individuals who enter Medicare at ages 20-62, and annual Medicare spending after first (usually partial) year of Medicare coverage
*   2262: main, but individuals who enter Medicare at ages 22-62
local sample 2262

* Sample definitions:
* - main (all ages, all years)
local sample_main 1
* - ages 22-62
local sample_2262 inrange(age_year_covstart_fill, 22, 62)
* - ages 51-52
local sample_5152 inrange(age_year_covstart_fill, 51, 52)
* - first full year
local sample_ffy (rfrnc_yr == init_yr + 1)

* Which controls to include in each regression
local controls_01 i.years_since_covstart
* local controls_02 i.years_since_covstart##i.fipscounty_init_g
* local controls_03 i.years_since_covstart##i.covstart_month
* local controls_04 i.years_since_covstart                      i.age_year_covstart_fill##i.male_init
* local controls_05 i.years_since_covstart##i.fipscounty_init_g i.age_year_covstart_fill##i.male_init
* local controls_06 i.years_since_covstart##i.covstart_month    i.age_year_covstart_fill##i.male_init
* local controls_07 i.years_since_covstart                      i.age##i.male_init
* local controls_08 i.years_since_covstart##i.fipscounty_init_g i.age##i.male_init
* local controls_09 i.years_since_covstart##i.covstart_month    i.age##i.male_init

* DI entrant age average Medicare spending and mortality
local clust county_mofd
foreach y in tot_pmt died_adj {
foreach x in age_year_covstart_fill {
  * local y tot_pmt
  * local x age_year_covstart_fill
  
  * Control variables
  local ctrl 01
  assert "`controls_`ctrl''" == "i.years_since_covstart"
  
  * Confirm sample
  assert "`sample'" == "2262"
  
  * Load estimates
  estimates use "$SSDIMed/results/estimates/x-i.`x'/x-i.`x'_y-`y'_controls-`ctrl'_cluster-`clust'_sample-`sample'.ster"
  
  * QC: confirm controls
  assert "`e(absvars)'" == "i.years_since_covstart"
  
  * Calculate y effect, CI for each level of x variable (_cons plus x-var effect)
  local addvar 0, 0, 0
  local names : colfullnames e(b)
  foreach name of local names {
    if "`name'" == "_cons" continue
    *local val = subinstr("`name'", "b.`=subinstr("`x'", "i.", "", .)'", "", .)
    *local val = subinstr("`val'",  ".`=subinstr("`x'", "i.", "", .)'", "", .)
    lincom _cons + `name'
    local addvar `addvar', `"`name' + _cons"', `r(estimate)', `r(se)'
  }
  regsave, addvar(`addvar') detail(all)
  
  * which variables to keep
  drop if var == "0"
  keep if strpos(var, "+ _cons") > 0
  
  * confidence intervals for added variables
  gen double t_ci = abs(invttail(`e(df_r)', (1 - 0.95) / 2))
  qui gen double ci_lower = coef - t_ci*stderr
  qui gen double ci_upper = coef + t_ci*stderr
  rename (coef ci_lower ci_upper) (`y' `y'_ci_lo `y'_ci_hi)
  keep var `y' `y'_ci_lo `y'_ci_hi clustvar
  
  * numeric level for the x-variable (a factor)
  gen `x' = substr(var, 1, strpos(var, "."))
  replace `x' = regexr(`x', "(b\.)|(\.)", "")
  destring `x', replace
  order `x', after(var)
  
  * Units
  * tot_pmt in $1,000
  * died in deaths per 10,000
  foreach var in `y' `y'_ci_lo `y'_ci_hi {
    if "`y'" == "tot_pmt"        replace `var' = `var'/1000
    if strpos("`y'", "died") > 0 replace `var' = `var'*10000
  }
  
  * Save
  tempfile `x'_`y'
  save ``x'_`y'', replace
  list
}
}

* DI incidence
local y entry
foreach x in age_year_covstart_fill {  
  * Load data
  use "$SSDIMed/data/analysis/county-month-age_entry_sample-main.dta", clear
  gisid fipscounty_firstnm covstart_month age_year_covstart_fill, missok
  
  * Confirm original sample size
  sum count
  assert r(sum) == 16240776
  qui sum age_year_covstart_fill
  assert r(min)==20 & r(max)==62
  
  * Select sample
  assert `sample_`sample'' == inrange(age_year_covstart_fill, 22, 62)
  keep if `sample_`sample''
  sum count
  assert r(sum) == 15790262

  * Collapse
  gcollapse (mean) nat_unemp_mean=unemp_rate_national_atapp  incidence=incidence_pop_age_atapp [aw = pop_age_atapp], by(`x') fast
  
  * Save
  tempfile `x'_`y'
  save ``x'_`y'', replace
  list, sep(0)
}



* ---------------------------------------------------------------------------------------
* Graph settings
* ---------------------------------------------------------------------------------------

* Set graph style to project default settings
*	Pass as arguments: height width (default 3.5in 6.5in)
do "$SSDIMed/scripts/_auxiliary/_project_grstyle.do" 3in 3in

local red "179 0 0"
local green "75 115 47"
local blue "48 84 150"
local bluegray "126 153 180"
local brown "140 045 004"
local orange "237 125 49"
local tan "210 180 140"

* Entry
local p 2
local color`p' `brown'
grstyle set color "`color`p''": p`p'
grstyle set symbol i: p`p'
grstyle set symbolsize 2.2pt: p`p'
grstyle set lpattern solid: p`p'
grstyle set linewidth 1.8pt: p`p'
* shaded region for confidence intervals
grstyle set color "`color`p''%20": p`p'area
grstyle set color "`color`p''%20": p`p'arealine
grstyle set linewidth vvthin: p`p'area

* Medical spending
local p 3
local color`p' black
grstyle set color "`color`p''": p`p'
grstyle set symbol i: p`p'
grstyle set symbolsize 2.2pt: p`p'
grstyle set lpattern solid: p`p'
grstyle set linewidth 1.8pt: p`p'
* shaded region for confidence intervals
grstyle set color "`color`p''%20": p`p'area
grstyle set color "`color`p''%20": p`p'arealine
grstyle set linewidth vvthin: p`p'area

* Mortality
local p 4
local color`p' `green'
grstyle set color "`color`p''": p`p'
grstyle set symbol i: p`p'
grstyle set symbolsize 2.2pt: p`p'
grstyle set lpattern solid: p`p'
grstyle set linewidth 1.8pt: p`p'
* shaded region for confidence intervals
grstyle set color "`color`p''%20": p`p'area
grstyle set color "`color`p''%20": p`p'arealine
grstyle set linewidth vvthin: p`p'area

* Arrows and lines
local p 5
grstyle set color gs8: p`p'
grstyle set lpattern solid: p`p'
grstyle set linewidth vvthin: p`p'


* -------------------------------------------------------------------------------------------------
* Figure: DI entrant characteristics, by age at Medicare entry
* -------------------------------------------------------------------------------------------------

* Load data
local x age_year_covstart_fill
use ``x'_tot_pmt', clear
merge 1:1 var using ``x'_died_adj', assert(match) nogen noreport
merge 1:1 `x' using ``x'_entry', assert(match using) nogen noreport
list, sep(0)
assert inrange(age_year_covstart_fill, 20, 62)

* xvar is age at entry, taken to be two years before Medicare coverage starts
cap drop xvar
gen double xvar = `x' - 2
assert inrange(floor(xvar), 22 - 2, 62 - 2)
qui glevelsof xvar
assert r(J) == (62 - 22 + 1)

* ------------------------------
* Panel (a): By age at Medicare entry: DI entrant cohort size
if 1 {

* How many major ticks?
* How much should range extend beyond ticks (as fraction of tick step size)
local y_ticks 6
local y_margin_lo 0
local y_margin_hi 0.3

* y1 tick lo, step size, format
local y1_tick_lo 0
local y1_step 200
local y1_fmt "%9.2gc"
local y1_unit ""

* Construct y labels and scales from parameters above
forvalues y = 1/1 {
  local y`y'_tick_hi  = `y`y'_tick_lo' + `y`y'_step' * (`y_ticks' - 1)
  local y`y'_range_lo = `y`y'_tick_lo' - `y`y'_step' * `y_margin_lo'
  local y`y'_range_hi = `y`y'_tick_hi' + `y`y'_step' * `y_margin_hi'
  di "y`y'_tick_hi:  `y`y'_tick_hi'"
  di "y`y'_range_lo: `y`y'_range_lo'"
  di "y`y'_range_hi: `y`y'_range_hi'"
  
  local y`y'_label ylabel(`y`y'_tick_lo'(`y`y'_step')`y`y'_tick_hi', axis(`y') format("`y`y'_fmt'") noticks nolabels labgap(1pt)) 
  local y`y'_scale yscale(`yalt' range(`y`y'_range_lo' `y`y'_range_hi') axis(`y') noline) 
  di `"y`y'_label:  `y`y'_label'"'
  di `"y`y'_scale:  `y`y'_scale'"'
}

* x-axis settings
local xlo 15.5
local xhi 62
local xscale xscale(range(`xlo' `xhi'))
local x_axis xlabel(20(5)60) `xscale' xtitle("Age at entry", margin(t=+.5))

* style and color for y-axis title/labels
local y1_title_style it
local y1_title_color gs8
local y1_label_style it
local y1_label_color gs8
local y1_label_size 9pt

local text_style sf
local text_color black

* y-axis labels placed above grid rule (build manually)
local y1_x `xlo'
local y1_place ne
local y1_just left

local y 1
local y`y'_la
foreach tick of numlist  `y`y'_tick_lo'(`y`y'_step')`y`y'_tick_hi' {
  if `tick' == `y`y'_tick_hi' {
    * white textbox to go behind label
    local y`y'_la `y`y'_la' text(`=`tick' + `y`y'_step'/100' `y`y'_x' " ", yaxis(`y') place(`y`y'_place') just(`y`y'_just') size(`y`y'_label_size') box lwidth(none) color(`y`y'_label_color') width(5) height(3) bcolor(white%70))
    *label
    local y`y'_la `y`y'_la' text(`tick' `y`y'_x' "{`y`y'_label_style':`=string(`tick', "`y`y'_fmt'")'`y`y'_unit'}", yaxis(`y') place(`y`y'_place') just(`y`y'_just') size(`y`y'_label_size') box lwidth(none) color(`y`y'_label_color') bcolor(none))
  }
  else {
    local y`y'_la `y`y'_la' text(`tick' `y`y'_x' "{`y`y'_label_style':`=string(`tick', "`y`y'_fmt'")'}", yaxis(`y') place(`y`y'_place') just(`y`y'_just') size(`y`y'_label_size') box lwidth(none) color(`y`y'_label_color') bcolor(none))
  }
}

* y-axis titles, horizontal orientation (build manually)
local y 1
local y1title text(`=`y`y'_tick_lo' + `y`y'_step'*(`y_ticks' - 1 + `y_margin_hi')' `xlo' "{`y`y'_title_style':Entrants per million residents}"
local y1title `y`y'title', yaxis(`y') place(`y`y'_place') just(`y`y'_just') box lwidth(none) color(`y`y'_title_color') bcolor(white) margin(b=+.5))
local y1title `y`y'title' ytitle(" ", axis(`y') margin(l=-8))

* y-axis settings
local y_axis1 `y1_label' `y1_la' `y1_scale' `y1title'

* Plot lines
local y1 incidence
local line_y1 (connected `y1' xvar, yaxis(1) pstyle(p2))

* Label line_y1
local text_1 text(310 42 "{`text_style':Monthly DI entry}"
local text_1 `text_1', yaxis(1) place(nw) just(right) box lwidth(none) width(30) color(`text_color') bcolor(white%90))
local xref = 44
qui sum `y1' if xvar == `xref'
local pci_1 (pci `=r(mean)' `=`xref'' `=r(mean) + 40' `=`xref' - 2', pstyle(p5))

* Eligibility changes
local text_grid text(`=`y1_tick_lo' + 0.62*`y1_step'' 52.5 "Eligibility" "relaxes", place(c) just(center) box lwidth(none) bcolor(white%94))
local xlines_grid xline(50) xline(55)

local legend legend(off)
local region graphregion(margin(t=5 b=-1))

twoway `line_y1' `pci_1', `x_axis' `y_axis1' `legend' `region' `text_1' `text_grid' `xlines_grid'

* Export graph
cap mkdir "$SSDIMed/results/figures"
graph export "$SSDIMed/results/figures/entry_age_cohort_count-`sample'.pdf", as(pdf) fontface("Source Sans Pro") fontdir("$SSDIMed/data/raw/fonts") replace

}


* ------------------------------
* Panel (b):  By age at Medicare entry: DI entrant Medicare spending (tot_pmt)
if 1 {

* How many major ticks?
* How much should range extend beyond ticks (as fraction of tick step size)
local y_ticks 6
local y_margin_lo 0
local y_margin_hi 0.3

* y1 tick lo, step size, format
local y1_tick_lo 8
local y1_step 2
local y1_fmt "%9.2gc"
local y1_sign "$"
local y1_unit " thousand"

* Construct y labels and scales from parameters above
forvalues y = 1/1 {
  local y`y'_tick_hi  = `y`y'_tick_lo' + `y`y'_step' * (`y_ticks' - 1)
  local y`y'_range_lo = `y`y'_tick_lo' - `y`y'_step' * `y_margin_lo'
  local y`y'_range_hi = `y`y'_tick_hi' + `y`y'_step' * `y_margin_hi'
  di "y`y'_tick_hi:  `y`y'_tick_hi'"
  di "y`y'_range_lo: `y`y'_range_lo'"
  di "y`y'_range_hi: `y`y'_range_hi'"
  
  local y`y'_label ylabel(`y`y'_tick_lo'(`y`y'_step')`y`y'_tick_hi', axis(`y') format("`y`y'_fmt'") noticks nolabels labgap(1pt)) 
  local y`y'_scale yscale(`yalt' range(`y`y'_range_lo' `y`y'_range_hi') axis(`y') noline) 
  di `"y`y'_label:  `y`y'_label'"'
  di `"y`y'_scale:  `y`y'_scale'"'
}

* x-axis settings
local xlo 15.5
local xhi 62
local xscale xscale(range(`xlo' `xhi'))
local x_axis xlabel(20(5)60) `xscale' xtitle("Age at entry", margin(t=+.5))

* style and color for y-axis title/labels
local y1_title_style it
local y1_title_color gs8
local y1_label_style it
local y1_label_color gs8
local y1_label_size 9pt

local text_style sf
local text_color black

* y-axis labels placed above grid rule (build manually)
local y1_x `xlo'
local y1_place ne
local y1_just left
local y1_width 25

local y 1
local y`y'_la
foreach tick of numlist  `y`y'_tick_lo'(`y`y'_step')`y`y'_tick_hi' {
  if `tick' == `y`y'_tick_hi' {
    * white textbox to go behind label
    local y`y'_la `y`y'_la' text(`=`tick' + `y`y'_step'/100' `y`y'_x' " ", yaxis(`y') place(`y`y'_place') just(`y`y'_just') size(`y`y'_label_size') box lwidth(none) color(`y`y'_label_color') width(`y`y'_width') height(4) bcolor(white%70))
    *label
    local y`y'_la `y`y'_la' text(`tick' `y`y'_x' "{`y`y'_label_style':`y`y'_sign'`=string(`tick', "`y`y'_fmt'")'`y`y'_unit'}", yaxis(`y') place(`y`y'_place') just(`y`y'_just') size(`y`y'_label_size') box lwidth(none) color(`y`y'_label_color') bcolor(none))
  }
  else {
    local y`y'_la `y`y'_la' text(`tick' `y`y'_x' "{`y`y'_label_style':`y`y'_sign'`=string(`tick', "`y`y'_fmt'")'}", yaxis(`y') place(`y`y'_place') just(`y`y'_just') size(`y`y'_label_size') box lwidth(none) color(`y`y'_label_color') bcolor(none))
  }
}

* y-axis titles, horizontal orientation (build manually)
local y 1
local y1title text(`=`y`y'_tick_lo' + `y`y'_step'*(`y_ticks' - 1 + `y_margin_hi')' `xlo' "{`y`y'_title_style':Medical spending}"
local y1title `y`y'title', yaxis(`y') place(`y`y'_place') just(`y`y'_just') box lwidth(none) color(`y`y'_title_color') bcolor(white) margin(b=+.5))
local y1title `y`y'title' ytitle(" ", axis(`y') margin(l=-8))

* y-axis settings
local y_axis1 `y1_label' `y1_la' `y1_scale' `y1title'

* Plot lines
local y1 tot_pmt
local xvar age_year_covstart_fill
local line_y1 (connected `y1' xvar, yaxis(1) pstyle(p3))
local ci_area (rarea `y1'_ci_lo `y1'_ci_hi xvar, sort pstyle(p3area))

* Label line_y1
local text_1 text(14.5 44 "{`text_style':Annual medical spending}"
local text_1 `text_1', yaxis(1) place(nw) just(right) box lwidth(none) width(50) color(`text_color') bcolor(white%90))
local xref = 46
qui sum `y1' if xvar == `xref'
local pci_1 (pci `=r(mean)' `=`xref'' `=r(mean) + 1' `=`xref' - 2', pstyle(p5))

* Eligibility changes
local text_grid text(`=`y1_tick_lo' + 0.62*`y1_step'' 52.5 "Eligibility" "relaxes", place(c) just(center) box lwidth(none) bcolor(white%94))
local xlines_grid xline(50) xline(55)

local legend legend(off)
local region graphregion(margin(t=5 b=-1))

twoway `ci_area' `line_y1' `pci_1', `x_axis' `y_axis1' `legend' `region' `text_1' `text_grid' `xlines_grid'

* Export graph
cap mkdir "$SSDIMed/results/figures"
graph export "$SSDIMed/results/figures/entry_age_cohort_tot_pmt_clust-`clust'-`sample'.pdf", as(pdf) fontface("Source Sans Pro") fontdir("$SSDIMed/data/raw/fonts") replace

}


* ------------------------------
* Panel (c):  By age at Medicare entry: DI entrant Medicare mortality (died)
if 1 {

* How many major ticks?
* How much should range extend beyond ticks (as fraction of tick step size)
local y_ticks 6
local y_margin_lo 0
local y_margin_hi 0.3

* y1 tick lo, step size, format
local y1_tick_lo 0
local y1_step 100
local y1_fmt "%9.2gc"
local y1_unit ""

* Construct y labels and scales from parameters above
forvalues y = 1/1 {
  local y`y'_tick_hi  = `y`y'_tick_lo' + `y`y'_step' * (`y_ticks' - 1)
  local y`y'_range_lo = `y`y'_tick_lo' - `y`y'_step' * `y_margin_lo'
  local y`y'_range_hi = `y`y'_tick_hi' + `y`y'_step' * `y_margin_hi'
  di "y`y'_tick_hi:  `y`y'_tick_hi'"
  di "y`y'_range_lo: `y`y'_range_lo'"
  di "y`y'_range_hi: `y`y'_range_hi'"
  
  local y`y'_label ylabel(`y`y'_tick_lo'(`y`y'_step')`y`y'_tick_hi', axis(`y') format("`y`y'_fmt'") noticks nolabels labgap(1pt)) 
  local y`y'_scale yscale(`yalt' range(`y`y'_range_lo' `y`y'_range_hi') axis(`y') noline) 
  di `"y`y'_label:  `y`y'_label'"'
  di `"y`y'_scale:  `y`y'_scale'"'
}

* x-axis settings
local xlo 15.5
local xhi 62
local xscale xscale(range(`xlo' `xhi'))
local x_axis xlabel(20(5)60) `xscale' xtitle("Age at entry", margin(t=+.5))

* style and color for y-axis title/labels
local y1_title_style it
local y1_title_color gs8
local y1_label_style it
local y1_label_color gs8
local y1_label_size 9pt

local text_style sf
local text_color black

* y-axis labels placed above grid rule (build manually)
local y1_x `xlo'
local y1_place ne
local y1_just left

local y 1
local y`y'_la
foreach tick of numlist  `y`y'_tick_lo'(`y`y'_step')`y`y'_tick_hi' {
  if `tick' == `y`y'_tick_hi' {
    * white textbox to go behind label
    local y`y'_la `y`y'_la' text(`=`tick' + `y`y'_step'/100' `y`y'_x' " ", yaxis(`y') place(`y`y'_place') just(`y`y'_just') size(`y`y'_label_size') box lwidth(none) color(`y`y'_label_color') width(5) height(3) bcolor(white%70))
    *label
    local y`y'_la `y`y'_la' text(`tick' `y`y'_x' "{`y`y'_label_style':`=string(`tick', "`y`y'_fmt'")'`y`y'_unit'}", yaxis(`y') place(`y`y'_place') just(`y`y'_just') size(`y`y'_label_size') box lwidth(none) color(`y`y'_label_color') bcolor(none))
  }
  else {
    local y`y'_la `y`y'_la' text(`tick' `y`y'_x' "{`y`y'_label_style':`=string(`tick', "`y`y'_fmt'")'}", yaxis(`y') place(`y`y'_place') just(`y`y'_just') size(`y`y'_label_size') box lwidth(none) color(`y`y'_label_color') bcolor(none))
  }
}

* y-axis titles, horizontal orientation (build manually)
local y 1
local y1title text(`=`y`y'_tick_lo' + `y`y'_step'*(`y_ticks' - 1 + `y_margin_hi')' `xlo' "{`y`y'_title_style':Deaths per 10,000}"
local y1title `y`y'title', yaxis(`y') place(`y`y'_place') just(`y`y'_just') box lwidth(none) color(`y`y'_title_color') bcolor(white) margin(b=+.5))
local y1title `y`y'title' ytitle(" ", axis(`y') margin(l=-8))

* y-axis settings
local y_axis1 `y1_label' `y1_la' `y1_scale' `y1title'

* Plot lines
local y1 died_adj
local line_y1 (connected `y1' xvar, yaxis(1) pstyle(p4))
local ci_area (rarea `y1'_ci_lo `y1'_ci_hi xvar, sort pstyle(p4area))

* Label line_y1
local text_1 text(265 42 "{`text_style':Annual mortality}"
local text_1 `text_1', yaxis(1) place(nw) just(right) box lwidth(none) width(30) color(`text_color') bcolor(white%90))
local xref = 44
qui sum `y1' if xvar == `xref'
local pci_1 (pci `=r(mean)' `=`xref'' `=r(mean) + 25' `=`xref' - 2', pstyle(p5))

* Eligibility changes
local text_grid text(`=`y1_tick_lo' + 0.62*`y1_step'' 52.5 "Eligibility" "relaxes", place(c) just(center) box lwidth(none) bcolor(white%94))
local xlines_grid xline(50) xline(55)

local legend legend(off)
local region graphregion(margin(t=5 b=-1))

twoway `ci_area' `line_y1' `pci_1', `x_axis' `y_axis1' `legend' `region' `text_1' `text_grid' `xlines_grid'

* Export graph
cap mkdir "$SSDIMed/results/figures"
graph export "$SSDIMed/results/figures/entry_age_cohort_died_clust-`clust'-`sample'.pdf", as(pdf) fontface("Source Sans Pro") fontdir("$SSDIMed/data/raw/fonts") replace

}







** EOF
